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ABSTRACT 

o 

I We simulate incompressible, MHD turbulence using a pseudo-spectral code. Our 

^ ■ major conclusions are as follows. 

1) MHD turbulence is most conveniently described in terms of counter propagating shear 
Alfven and slow waves. Shear Alfven waves control the cascade dynamics. Slow waves 
play a passive role and adopt the spectrum set by the shear Alfven waves. Cascades 
\ composed entirely of shear Alfven waves do not generate a significant measure of slow 

waves. 

| 2) MHD turbulence is anisotropic with energy cascading more rapidly along k± than 

along fen, where k± and fey refer to wavevector components perpendicular and parallel 
to the local magnetic field. Anisotropy increases with increasing k±_ such that excited 
modes are confined inside a cone bounded by feii oc k\_ where 7 < 1. The opening angle 
of the cone, 0(fej_) oc k^ 1 , defines the scale dependent anisotropy. 

3) The ID inertial range energy spectrum is well fit by a power law, E(k±) oc with 
a > 1. 

4) MHD turbulence is generically strong in the sense that the waves which comprise it 



o 



suffer order unity distortions on timescales comparable to their periods. Nevertheless, 
turbulent fluctuations are small deep inside the inertial range. Their energy density is 
less than that of the background field by a factor 0( q ~ 1 )/( 1 ~t) <^ \ 
5) MHD cascades are best understood geometrically. Wave packets suffer distortions 



as they move along magnetic field lines perturbed by counter propagating waves. Field 
lines perturbed by unidirectional waves map planes perpendicular to the local field into 
each other. Shear Alfven waves are responsible for the mapping's shear and slow waves 
for its dilatation. The amplitude of the former exceeds that of the latter by 1/G(fej_) 
which accounts for dominance of the shear Alfven waves in controlling the cascade 
dynamics. 

6) Passive scalars mixed by MHD turbulence adopt the same power spectrum as the 
velocity and magnetic field perturbations. 

7) Decaying MHD turbulence is unstable to an increase of the imbalance between the 
flux of waves propagating in opposite directions along the magnetic field. Forced MHD 
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turbulence displays order unity fluctuations with respect to the balanced state if excited 
at low k± by 5(t) correlated forcing. It appears to be statistically stable to the unlimited 
growth of imbalance. 

8) Gradients of the dynamic variables are focused into sheets aligned with the magnetic 
field whose thickness is comparable to the dissipation scale. Sheets formed by oppositely 
directed waves are uncorrelated. We suspect that these are vortex sheets which the mean 
magnetic field prevents from rolling up. 

9) Items (l)-(6) lend support to the model of strong MHD turbulence put forth by 
Goldreich & Sridhar (GS). Results from our simulations are also consistent with the 
GS prediction 7 = 2/3, as are those obtained previously by Cho & Vishniac. The sole 
notable discrepancy is that ID energy spectra determined from our simulations exhibit 
a 3/2, whereas the GS model predicts a = 5/3. Further investigation is needed to 
resolve this issue. 



1. INTRODUCTION 

Most of the baryonic matter in the universe has such high electrical conductivity that magnetic 
fields diffuse very slowly through it. Thus fluid motions and motions of magnetic field lines are 
closely coupled. Large scale motions are generally turbulent, and incompressible MHD is the 
simplest approximation under which these complex coupled motions can be investigated. 

The inertial range of MHD turbulence is an essential ingredient in a variety of astronomical 
phenomena. Cosmic rays are scattered by inertial range magnetic field fluctuations. This affects 
both their propagation and their acceleration in shock fronts (Blandford & Eichler 1987; Berezinskii 
1990; Chandran 2000). Reconnection of magnetic field lines is an important ingredient of flare 
activity and dynamo action. The rate at which it proceeds seems likely to depend upon the small 
scale structure of magnetic field lines (Parker 1979; Lazarian & Vishniac 1999). The scintillation 
of small angular diameter radio sources due to scattering by electron density fluctuations is almost 
certainly related to inertial range MHD turbulence (Higdon 1984; Rickett 1990). 

The organization of this paper is as follows. Relevant properties of MHD waves are described 
in §2. In §3 we introduce selected analytical models for the inertial range of MHD turbulence. The 
strategy we follow in designing our simulations is set forth in §4. Results from these simulations 
are presented in §5. In §6 we interpret our results and compare them with results from prior 
investigations. Technical details of our simulation method are relegated to the appendix. Two 
subsections merit special mention. 

A comparison between our simulations and those of earlier workers is given in §6.5. Here we 
merely point out that our simulations are most closely related to those presented by Cho & Vishniac 
(2000). These authors demonstrate that under nearly isotropic forcing, MHD turbulence develops 
a scale dependent anisotropy which increases with increasing wave number in the manner suggested 
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by Goldreich & Sridhar (1995) and Goldreich Sz Sridhar (1997). Henceforth we refer to Goldreich 
& Sridhar (1995) and Goldreich & Sridhar (1997) separately as GSI and GSII, and together as GS. 
Our simulations differ from those of Cho & Vishniac (2000) in that they are excited anisotropically 
so that we can study the deep inertial range of MHD turbulence. 

Our most perplexing result, the shallow slope we find for the ID energy spectrum, is discussed in 
§6.6. Unfortunately we cannot offer a definitive explanation. This will require further investigation. 



2. BACKGROUND MATERIAL 
2.1. Basic Equations 

The equations which govern magnetohydrodynamics (MHD) are written below using notation 
defined in Table 1. 

(B 2 \ 1 
P + ^) + 4^ B ' VB + ^V 2 v, (1) 

d t B = V x (v x B) + z/ B V 2 B, (2) 

d t p + V • (pv) = 0, (3) 

V • B = 0. (4) 
The concentration of a passive scalar advected by the fluid evolves according to 

d t c + V- (cv) = z/ c V 2 c. (5) 

We simplify equations (l)-(5) for applications in this paper. 3 Incompressibility is assumed 
throughout, so we set p = 1 and define the total pressure P = p + B 2 /8tt. The magnetic field 



3 Further steps are taken in §A.2 to cast these equations in a form suitable for computation. 



Table 1. Notation 



v 


fluid velocity 


B 


magnetic field 


p 


fluid density 


V 


fluid pressure 


c 


passive scalar concentration 


v v 


momentum diffusivity 




magnetic diffusivity 




passive scalar diffusivity 
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is measured in velocity units by b = B/\/47r. Each diffusive term is replaced by a n'th order 

hyperdiffusivity with the same coefficient v n . With these modifications, equations (l)-(5) transform 
to 

8 t v = -v ■ Vv - VP + b ■ Vb + u n V 2n v, (6) 

8 t b = -v ■ Vb + b ■ Vv + z/ n V 2n b, (7) 

V • v = 0, (8) 

V • b = 0. (9) 

d t c + vVc = v n V 2n c. (10) 

To relate P to v and b, we take the divergence of equation (6) which yields 

V 2 P = Vb : Vb - Vv : Vv. (11) 

Thus 



f d 6 x' Vv : Vv - Vb : Vb , , 

p = / n—- n 1 ■ 12 

J Air |x' - x| 

2.2. Regimes 

We decompose the magnetic field into a uniform part plus fluctuations; 

b = (b) + Ab. (13) 

The Alfven speed is defined by t^z = (b); va is taken to be constant in space and in time as is 
consistent with flux conservation. The energy densities of the mean magnetic field, the velocity 
field, and the magnetic fluctuations are denoted by E^, E v , and £^Ab- The parameter 

M= p , (14) 

which measures the relative importance of the fluctuations compared to the uniform field, deter- 
mines the character of MHD turbulence. 

MHD turbulence with small jjl can be described in terms of interacting waves. Kinetic and 
potential energy are freely interchanged so E v and Ej\b have comparable magnitudes. Wavemode 
turbulence is the principal subject of this thesis. Analytic scalings are presented in §3 to provide 
an intuitive feel for its dynamics. Results from our simulations are described in §5 and discussed 
in S6. 
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2.3. Linear Waves In Incompressible MHD 

Linear perturbations about a uniform background magnetic field can be decomposed into shear 
Alfven and pseudo Alfven waves. The pseudo Alfven wave is the incompressible limit of the slow 
magnetosonic wave. 4 As is well known, both waves conform to the dispersion relation 

u? = v\k 2 z . (15) 

Eigenvectors for these modes take the form 

v^(k, t) = a(k) expik • (x =F w^tz), b^k, t) = Ta(k) exp ik • (x =F VAtz), (16) 

vs(k, t) = s(k) expik • (x =F VAtz), bs(k, t) = =Fs(k) exp ik • (x =F VAtz), (17) 
where the unit polarization vectors are defined by 

k x z z — (k • z)k , . 

a= ; -, s= , '—. (18) 

[1 - (k-z) 2 ] 1 ^' [1- (k-z) 2 ]V2 

We note that k, s, and a form a right-hand triad. 

MHD turbulence is anisotropic with power cascading more rapidly to high k± than to high k z . 
In the limit k± » k z , s —>■ z; displacements associated with slow modes align along the unperturbed 
magnetic field. 



2.4. Elsasser Variables 

The Elsasser transformation 

Wf = vaz + v — b w | = —vaz + v + b (19) 

applied to equations (6) and (7) with v n = brings out the two wave characteristics 

9tW| + VAd z w^ = — Wj • Vw| — VP, (20) 
<9 t W| — vaQz^i = — Wf • Vwj_ — VP, (21) 

where from equation (12), 

d 3 x' Vw| : Vwj 



/ 



4tt 



(22) 



Linear waves propagate at the Alfven speed va either parallel (wj) or anti-parallel (wjj to the 
direction of the background magnetic field. 



4 In the limit of incompressibility the fast magnetosonic wave has infinite phase velocity and cannot be excited. 
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2.5. Collisions Between Wave Packets 

A disturbance of the background field may be decomposed into upward (w-f) and downward 
(w|) propagating wave packets. In the special case of unidirectional propagation either wj = 
or wj_ = 0, and an arbitrary nonlinear wave packet is an exact solution of the equations of 
incompressible MHD (Parker 1979). To prove this, take the divergence of the equation for the 
nonzero w. This yields V 2 P = which, since it applies globally, implies VP = 0, and hence 
that the wave packet propagates without distortion. 5 An important corollary is that nonlinear 
distortions occur only during collisions between oppositely directed wave packets. 

Collisions are constrained by the conservation laws of energy, 

E = \j d 3 x(|v| 2 + |b| 2 ), (23) 

and cross helicity, 

/= I J d 3 xvh. (24) 

These conservation laws follow directly from equations (6)-(9) in the limit that v n = 0. As a 
consequence, energy is not exchanged between colliding wave packets. A short proof follows. 

Take the dot product of equations (20) and (21) with wj and wj^, respectively. The advective 
and pressure gradient terms reduce to total divergences. This establishes that 

!/*x|w t |»-£/*,|w 1 |'-0 (25) 

provided wj and wj, either vanish at infinity or satisfy periodic boundary conditions. /,From the 
defining equations for the Elsasser variables, we obtain |wf| 2 /2 = |v| 2 + |b| 2 for wj_ = and 
| Wi | 2 /2 = |v| 2 + |b| 2 for w T = 0. Thus 

Er = \J d 3 *|w T | 2 and E l = ^J d 3 x\ Wl \ 2 (26) 

are the energies of isolated upward and downward propagating wave packets. This completes the 
proof that wave packet collisions are elastic. 



2.6. Wave Packets Move Along Field Lines 

To lowest nonlinear order in the wave amplitudes, distortions suffered in collisions between 
oppositely directed wave packets arise because each packet moves along field lines perturbed by the 



5 This conclusion remains valid for our simulations which are carried out in a computational box and employ 
periodic boundary conditions. 
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other. The proof follows directly from equation (20) written to second order in the amplitudes of 
the Wj and fields. With the aid of equation (22), it can be shown that 



D^ 



(i) 



d 3 x > V£« : Vw} 



47T 



X' 



X 



Here 



and 



5_ _9_ 

<9t <9z 



x(x ,t) = x + £(x ,t). 



(27) 



(28) 



(29) 



The Lagrangian displacement, £, connects the Lagrangian coordinate of a fluid particle, xo, to its 
Eulerian coordinate, x. 



Several steps are needed to establish equation (27). In terms of £, 



dt 



V A Z + V A 



XO 



dz 



(30) 



To first order in the amplitudes of wj and wj_, we may replace x by x in the definition of £ and 
write 



(i) 



at 



a) 



dz 



It then follows from equations (19) and (31) that 



w 



(i) 



(31) 



(32) 



The final step is to verify that the linear operator passes through the integral sign and changes 



£^ to wr ; while leaving the rest of the integrand unaltered. 



i 



w 



(2) 



Equation (27) has a simple interpretation. Consider an upward moving wave packet for which 
prior to its interaction with downward moving waves. Subsequent to this interaction 



suppose that ^ at fixed Zj = z — v^t is changed by A^ 1 ). Then as a function of z^ 

tfV VA£« : Vwf 1 



Awf^-A^-Vw^-V J 



4tt 



x' 



X 



(33) 



The first term on the right-hand side of this equation is the perturbation that would result from 
the unconstrained displacement of wj 1 ^ by A£^ . The second term constrains the perturbation to 

preserve V ■ w| =0. Since magnetic field lines are frozen in the fluid, we conclude that, at least 
to second order, wave packets follow magnetic field lines. 

Two points are worth stressing in connection with equation 33. 



6 This term arises from the gradient of the total pressure. 
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• Downward propagating waves contribute the entire A£ ( ' since it is measured at fixed = 
z — VAt- This is consistent with the general rule that only oppositely directed wave packets 
interact. 

• The turbulent energy cascade is associated with the shear of the field. Uniform dis- 
placements, which arise as a consequence of the sweeping of small disturbances by larger ones, 
do not contribute to the transfer of energy across scales. 



The proof in this section has been couched in Eulerian coordinates. A technically simpler 
version in Lagrangian coordinates is given by Sridhar k, Goldreich (1994). It consists of demon- 
strating that the third order Lagrangian density for incompressible MHD vanishes when written in 
terms of the transverse components of the displacement vector. Although simpler technically, the 
Lagrangian based result is more subtle conceptually. Its proper interpretation is provided in GSII. 



3. MHD CASCADES 

A variety of models have been proposed for MHD turbulence. They share the common feature 
that energy cascades from lower to higher wave number. 



3.1. The Iroshnikov-Kraichnan Model 

The standard model is that due to Iroshnikov (1963) and Kraichnan (1965). Kraichnan's 
derivation of the IK spectrum relies on the fact that only oppositely directed waves interact in 
incompressible MHD. It assumes explicitly that the turbulence is isotropic and implicitly that the 
dominant interactions are those which couple three waves. 

The above assumptions imply that the cascade time across scale A is 

fc-f 24 )'-- (34) 
Setting v\jt c equal to the dissipation rate per unit mass, e, then yields 

v\ ~ (ev A X) 1/4 , (35) 
which corresponds to the ID power spectrum 7 



7 Because the IK cascade is isotropic, it is sufficient to specify its ID power spectrum. 
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Nonlinearity is measured by % ~ (v\/va), where N ~ x is the number of wave periods in t c ; 



Since x decreases with decreasing A, only dissipation limits the length of the IK inertial range. 

The IK model is flawed because the assumption of isotropy is inconsistent with the frequency 
and wavevector closure relations that resonant triads must satisfy (Shebalin, Matthaeus, & Mont- 
gomery 1983). These take the form 

UJl + Ll> 2 = W3) (38) 

k!+k 2 = k 3 . (39) 

But since u = va \ k z | , equation (38) and the z component of equation (39) yield the set 

\ki z \ + \k 2z \ = \k 3z \ (40) 
ku + k 2z = k 3z . (41) 

Because nonlinear interactions can only occur between oppositely directed waves, the 3-mode cou- 
pling coefficient vanishes unless waves 1 and 2 propagate in opposite directions. In that case, 
equations (40) and (41) imply that either k\ z or k 2z must vanish. Since one of the incoming waves 
has zero frequency, 3- wave interactions do not cascade energy along k z . 



3.2. Intermediate MHD Turbulence 

GSII propose an anisotropic MHD cascade based on scalings obtained from 3- wave interactions. 
It represents a new form of turbulence, which they term intermediate, because it shares some of 
the properties of both weak and strong turbulence. Although individual wave packets suffer small 
distortions in single collisions, interactions of all orders make comparable contributions to the 
perpendicular cascade. 8 

To derive the scaling relations for the intermediate cascade, we repeat the steps carried out 
in §3.1 for the IK model, but with Aj_ in place of A and Ay held constant. Here A^ and Ay are 
correlation lengths in directions perpendicular and parallel to the local magnetic field. Thus 

^~(^) 2 -- < 42 > 



Setting e ~ v\ /t c , we find 



^~(^) 1/4 , (43) 



8 This is a controversial claim. 
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and 

E{k±) Jj^l. (44) 

Besides being anisotropic, the intermediate MHD cascade differs from the IK cascade in another 
important respect. The strength of nonlinear interactions, as measured by 

increases along the cascade. Thus, even in the absence of dissipation, the intermediate cascade has 
a finite inertial range. This suggests that a strong form of MHD turbulence must be the relevant 
one for most applications in nature. 



3.3. Strong MHD Turbulence 

A cascade for strong MHD turbulence is described in GSI. Its defining property is that MHD 
waves suffer order unity distortions on time scales comparable to their periods. This state is referred 
to as one of critical balance. Motivation for the hypothesis of critical balance is given in Goldreich 
& Sridhar (1995, 1997) and summarized below. Our discussion of intermediate turbulence shows 
that x increases if it is less than unity. However, it cannot rise above unity since the frequency 
spread of the wave packets which emerge following a strong collision must satisfy the frequency- 
time uncertainty relationship. Gruzinov (2000) provides a more physical explanation for the upper 
bound on \- He points out that for x S> 1, 2D motions of scale A^ in planes perpendicular to 
the local magnetic field are uncoupled over separations greater than \\/x along the field direction. 
Thus, during a time interval of order Aj_/da ± ~ \\/vaX these motions reduce x to order unity. 

Crtical balance, and the assumption of a constant energy flux along the cascade as expressed 

by 

e~-*i (46) 
Aj_ 

imply 

Ay oc X 2 / 3 . (47) 

Although there is a parallel cascade of energy in strong MHD turbulence, the degree of anisotropy 
increases along the cascade. 

Let us assume vl ~ va and isotropy on scale outer scale L. Then the 3D energy spectrum of 
strong MHD turbulence takes the form 

E <^>~z4^(lf )' (48> 
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where f(u) is a positive symmetric function of u with the properties that f(u) ps 1 for |u| < 1 and 
f{u) negligibly small for \u\ S> 1. The power spectrum is flat as a function of k\\ for fen < fc^ 3 L -1 / 3 
because the velocity and magnetic perturbations on transverse scale A?^ 1 arise from independent 
wave packets whose lengths Am ~ Af LVs. The ID perpendicular power spectrum obtained from 
equation (48) reads 

E{kl) ~ l4f' <49> 

Thus the spectrum of strong MHD turbulence is an anisotropic version of the Kolmogorov (1941) 
spectrum for hydrodynamic turbulence. 

Inertial range velocity differences and magnetic perturbations across perpendicular scale Aj_ 
satisfy 

/A ± \ 1/3 

v\ ± ~ b x± ~ I — J v A . (50) 



Thus even though the turbulence is properly classified as strong, deep in the inertial range magnetic 
field lines are nearly parallel across perpendicular separations Aj_ and nearly straight along parallel 
separations An; differential bending angles are of order (A^/L) 1 / 3 ~ (An/L) 1 / 2 . 



3.3.1. Parallel Cascade 

It is interesting to examine the frequency changing interactions that drive the parallel cascade. 
Referring back to the intermediate cascade, we know that 3-wave interactions do not change fre- 
quencies. However, interactions involving more than 3-waves can. For example, frequency changes 
arise in 4-wave interactions of the form 

OJl + Ll> 2 + W 3 = W4, (51) 

ki+k 2 + k 3 = k 4 , (52) 

where k\ z and have the same sign and W3 = fA|^3z| = (Ng & Bhattacharjee 1996, GKII). The 
parallel cascade they give rise to proceeds at a rate which is smaller than that of the perpendicular 
cascade by a factor of order \- Because strong MHD turbulence is characterized by x ~ 1> it has 
a significant parallel cascade. 



3.3.2. Field Line Geometry 

MHD turbulence is best understood geometrically. Field lines perturbed by waves propagating 
in one direction define two-dimensional mappings between xy planes separated by distance z. Shear 
Alfven waves dominate the shear and slow waves the dilatation of these mappings. The magnitude 
of the shear exceeds that of the dilatation by a factor of order Am/Aj_ ~ (-^/AjJ 1 / 3 3> 1. These 
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mappings describe the distortion that counter propagating waves would suffer if they moved at 
uniform speed along the perturbed field lines. The dominance of the shear over the dilatation 
explains why shear Alfven waves control the perpendicular cascades of both types of wave. 

The recognition that MHD waves tend to follow field lines is essential to understanding their 
turbulent cascades. Figure 1 provides a visual illustration of how this works. The left-hand panel 
displays a snapshot of field lines perturbed by downward propagating waves. In the right-hand 
panel we follow the evolution of a horizontal pattern as it propagates from the bottom to the top 
following these lines. The distortion of the initially circular bullseye is principally due to the shear 
in the two-dimensional mapping defined by the perturbed field lines. The cascade time on the scale 
of the initial pattern is that over which the shear grows to order unity. 

This geometrical picture requires two qualifications. The first is that the propagation speed of 
MHD waves is not exactly constant but varies with the strength of the local magnetic field. Pressure 
perturbations associated with slow waves are balanced by perturbations of magnetic pressure. The 
resulting perturbations in propagation speed, of order v\ ± , contribute to the nonlinear cascade. 
Over one wave period they lead to fractional distortions of order v\ ± /va ~ Aj_/A| <C 1. Thus they 
are properly ignored. The second qualification is that MHD waves do not exactly follow field lines. 
The extent to which this effects their cascade remains to be quantified. 

The parallel cascade may also be viewed in geometrical terms. Consider an upward propagating 
wave packet of length Ay and width Aj_ which is being distorted by downward moving wave packets 
of similar scale. Correlations along the parallel direction are shortened because the front and back 
of the wave packet undergo different 2-dimensional mappings. This happens because the upward 
propagating packet distorts each downward going packet as it passes through it. This distortion is 
of order \. For strong MHD turbulence \ ~ 1 which accounts for its significant parallel cascade. 

Incidentally, the geometrical picture also aids the interpretation of results from perturbation 
theory. For example, the 3-wave resonant interactions which dominate the perpendicular cascade 
and the 4-wave resonant interactions which cause the lowest order frequency changes each depend 
upon the amplitudes of modes with k z = 0. This is because the shear in the mapping between xy 
planes separated by Az is proportional to the displacement amplitudes of modes with k z < 1/Az. 
Perturbation theory corresponds to the limit of vanishing cascade strength in which shears of order 
unity are achieved in the limit of infinite separation along the z-axis. 

3.3.3. Relation to 2D Hydrodynamic Turbulence 

Magnetic field lines possess a tension which makes them ill-disposed to bend, but they are easily 
shuffled. This accounts for the 2D character of MHD turbulence. It also prompts an inquiry about 
the relation of MHD turbulence to 2D hydrodynamic turbulence. Each fluid elements conserves 
its vorticity in inviscid 2D hydrodynamics. This results in a direct cascade of enstrophy (vorticity 
squared) toward high k± and an inverse cascade of energy toward small k± (Lesieur 1990). As we 




Fig. 1. — Wavepacket Distortion Through Fieldline Wander. The left-hand panel displays a sample 
of field lines perturbed by downward propagating waves. The distortion of an originally circular 
bullseye pattern as it moves upward following these field lines is shown in the right-hand panel. 
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now demonstrate, MHD turbulence does not share these characteristics. 

The vorticity equation in MHD, obtained by taking the curl of equation (6) with u n = 0, reads 

° (V ^ X V) = V x [v x (V x v) - b x (V x b)] . (53) 

We concentrate on shear Alfven waves since they dominate the field aligned vorticity for nearly 
perpendicular cascades. Scaling the terms in equation (53) shows that 



dt \ ± 



(54) 



Thus changes on the cascade time scale; it is not even approximately conserved. Consequently, 
there is no enstrophy constraint to prevent energy from cascading toward larger k± . 



4. SIMULATION STRATEGY 

What follows is a comprehensive discussion of the techniques used in our simulations. Technical 
aspects of the spectral method are presented in the Appendix. 

4.1. Spectral Wave Mode Decomposition 

Separation of v(k) and b(k) into upward and downward propagating components is accom- 
plished by forming Fourier coefficients of the Elsasser variables w-f(k) and Wj(k) according to 

w T (k) = v(k) - b(k) w 4 (k) = v(k) + b(k). (55) 

Projections of w-j-(k) and w^(k) along the polarization directions of the linear incompressible MHD 
eigenmodes given by equation (18) yield amplitudes of upward and downward propagating Alfven 
and slow waves. In obvious notation 

A T (k) = a-w T (k) = iJ A (k)-Mk) 5 T (k)^s-W T (k) = %(k)-6 5 (k) (56) 

A t (k) = a • w^k) = i; A (k) + Mk) S x (k) = s ■ Wi(k) = 5 s (k) + b s (k) (57) 

where 

u A (k) = a- v(k) 6 A (k) = a-b(k) u s (k) = s • v(k) b s (k) = s • b(k). (58) 

The eigenmode frame, (k, s, a), is tied to the direction of the mean field, bo = z, to which the 
local field direction, b, is inclined by an angle 6 ~ vl ± /va- Consequently, our method for spectral 
decomposition erroneously mixes Alfven and slow modes. However, for nearly transverse cascades 
the mixing is only of order ^ 2 < 1. 
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Field line tilt also causes k r and k z to differ from fej_ and k\\ ; 

k r = cos 6k± + sin #fey k z = — sm Ok _\_ + cos 6k\\. (59) 

Thus k r ps fcj_[l + 0(# 2 )]. However, fc z ps fey + #&:_|_ ps #A;_|_, where the final relation applies because 
the degree of anisotropy increases with increasing fej_ along MHD cascades. Thus, k± can be 
represented to acceptable accuracy by k r . However, fey cannot be obtained from k z . Henceforth, we 
treat as equivalent k r and fej_ and L x = L y and L±. However, we are always careful to distinguish 
k z from fey and to note that 

k z ps ^k ± . (60) 

va 



4.2. Power Spectra 

Three-dimensional power spectra of field quantities, £30 (k), are azimuthally symmetric func- 
tions of k r = fe^x + k y y at fixed k z . 9 Accordingly, we define the 2D integrated power spectrum 
by 

E 2D (k r , k z ) = k r d(j)E 3D (k). (61) 

J 



It is important to note that Ez£>(k r , k z ) is not equivalent to E2T>(k±, fey). Moreover, the 
latter cannot be derived from the former. This shortcoming is due to the failure of the spectral 
decomposition procedure described in §4.1 to determine fey. It means that the 2D power spectrum 
is not a useful quantity. 10 However, we make so much use of the ID integrated power spectrum 
defined by 

/oo 
dk z E 2D (k r , k z ), (62) 
-00 

that henceforth we drop the subscript ID. 



4.3. Structure Functions 

The 3D behavior of MHD turbulence is best captured in real space using second order structure 
functions tied to the local magnetic field. We define transverse and longitudinal structure functions 
for the vector field U by 

SFTu(x±) =< [U(x' + x ± ) - U(x')] • [U(x' + x ± ) - U(x')] >, (63) 

9 In any specific realization this is true only in a statistical sense. 
10 We obtain 2D information from structure functions. 
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where xj_ • b = 0, and 



with f(xy) = J c 




length from x'. Averaging over x' is done with random volume sampling. Since the vector fields of 
interest possess statistical axial symmetry about b, we include an axial averaging of the direction 
of xj_ at fixed x± = |xj_| in the computation of SFTu(x±). 



The anisotropy of MHD turbulence complicates the discussion of constraints on the timestep 
and hyperviscosity. Accordingly, we begin by discussing the simpler case of spectral simulation of 
isotropic hydrodynamic turbulence. These constraints are summarized in Figure 2. 



We assume the Kolmogorov scaling. Given velocity j)i on outer scale L, inertial range velocity 
differences across A < L scale as v\ ~ (A/L) 1 / 3 ^ down to inner scale £ ~ (z^i/i^L 2 "^ 1 ) 3 ^ 6 "^ 2 ^. 

Four conditions constrain the values of the timestep, At, and hyperviscosity, v n , suitable for 
a spectral simulation of isotropic hydrodynamic turbulence. Each refers to the behavior of modes 
with the largest wavevectors, ku- We express these constraints in terms of the dimensionless 
variables At = VLkMAt and V = {y n k 2 ^) / {vLku)- 

• Conditions 1) and 2) are concerned with computational accuracy. 

1. Advection by outer scale eddies gives rise to fractional changes of order v^kuAt in the 
Fourier components of the smallest scale modes during one timestep. 11 Accurate com- 
putation requires 



which is the spectral equivalent of the Courant condition in real space. 

2. Hyperviscosity causes a fractional decay of order v n kj^At in the amplitudes of the smallest 
scale modes during a single timestep. Thus 



4.4. Timestep And Hyperviscosity 



4-4- 1- Isotropic Hydrodynamic Turbulence 



At < 1 



(65) 




(66) 



ii 



Changes caused by interactions which are local in Fourier space are smaller by a factor {UmL) 



-1/3 
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• Conditions 3) and 4) are required to maintain stability. 

3. This constraint depends upon the algorithm used to advance the variables in time. Inte- 

gration with RK2 results in an unphysical transfer of energy from large to small scale 
modes. Consider ID uniform advection at speed vl of the single Fourier mode v{k,M,t). 
RK2 yields v(k M ,At) = [1 - iv L k M At - (v L k M At) 2 /2]v(k M ,0). Thus \v(k M , At)\ 2 = 
[1 + (vLkMAt) 4 /A]\v(kM , 0)| 2 . In order that hyperviscosity maintain stability, 

At<{AV n f\ (67) 

4. A turbulent cascade transfers energy from large to small scales where it is dissipated by 

viscosity. Spectral simulations of turbulence must include a mechanism which is able 
to dispose of the energy carried by the cascade before it reaches ku- Otherwise it 
would reflect back to smaller k and the high k Fourier modes would approach energy 
equipartition with those of lower k. 12 Hyperviscosity suffices provided the inner scale it 
sets is larger than the grid resolution. This requires 

V>{k M L)- 1/Z = {irN/2)- l l\ (68) 

Dealiasing also involves a loss of energy and can stabilize simulations run with a suf- 
ficiently small timestep even in the absence of hyperviscosity. Further investigation is 
needed to clarify the manner in which energy is lost due to dealiasing. 



4-4-2. Anisotropic Magnetohydrodynamic Turbulence 

We restrict our discussion of MHD turbulence to cases in which the energy in the mean 
magnetic field greatly exceeds that in the kinetic and magnetic fluctuations. As discussed in Chapter 
3, analytic arguments indicate that the Kolmogorov scaling is obeyed in planes perpendicular 
to the mean magnetic field. Thus constraints on the timestep and hyperviscosity deduced in 
§4.4.1 for hydrodynamic turbulence pertain to MHD turbulence in the xy plane provided we take 
At = v L± k M± At and V = (v n k 2 %J/(v L± k M± )- 

Different constraints arise from motion along the direction of the mean magnetic field. Strong 
MHD turbulence is anisotropic with energy cascading more rapidly along k± than along fey . Analytic 
arguments imply that the anisotropy at perpendicular scale kj 1 is determined by the condition that 
the nonlinearity parameter \ = ( v \k±)/(vAk\\) ~ 1, where ku is the wavevector component in the 
direction of the local magnetic field. It is important to maintain the distinction between k\\ and k z . 
As discussed in §4.1, k z m (vL ± /vA)k± ~ {k±L±) 1 / 3 k» . 



The energy per Fourier mode scales as k ' in Kolmogorov turbulence. 
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Fig. 2. — Timestep and Hyperviscosity. We illustrate constraints on the dimensionless timestep 
and hyperviscosity as described in §4-4- The constraint given by equation (68) is not shown because 
it depends upon an additional parameter. Allowed choices lie in the unshaded part of the figure. 
Each plotted point represents values used in an individual simulation. 
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We control the value of (vl ± L z )/(vaL±) in our simulations. Typically this quantity is set 
somewhat larger than unity in order to ensure that the largest scale structures cascade on a time 
scale shorter than the Alfven crossing time L z /va- As a consequence, 

VAkM z ~ VaL ^~ v L± k M± < v L± k M± . (69) 

VL ± ^z 

After this preparation, we are ready to examine the constraints placed on At and v n by evolution 
in the z direction. 

1. Advection in the z direction is dominated by propagation at the Alfven speed since vl z <C v a 

in our simulations. Thus computational accuracy demands tvi&M z At < 1. Since t>A^M z ^ 
VL ± kM ± , this constraint is less severe than that imposed by equation (65). 

2. Hyperdiffusivity is not important in the z direction because ku z <C kM ± and we use a scalar 

hyperdiffusivity. 

3. Integration with RK2 leads to an unphysical transfer of energy from large to small scales due 

to advection at the Alfven speed. Provided equation (67) is satisfied, this does not cause any 
difficulty because VAku z ^5 fL ± &M ± - 

4. The maximum wave number in the z direction, ku z , must be larger than that at the inner scale 

of the cascade. From equation (60) 



As mentioned above, the factor preceding kM ± L± is typically larger than unity. Thus adequate 
resolution along z generally requires N z > N±. 



We carry out simulations of both forced and decaying MHD turbulence. Amplitudes of Fourier 
modes within 3 lattice units of the origin, normalized wave vector |s| < 3, are incremented at each 
timestep in simulations of forced turbulence and assigned initial values in simulations of decaying 
turbulence. Each component of these amplitudes receives an addition of a complex number with 
random phase and absolute value drawn from a Boltzmann distribution with specified mean subject 
to the constraint that k- v(k) = and k-b(k) = 0. Thus we are forcing both velocity and magnetic 
fluctuations. Forcing v alone would artificially correlate the power received by Wj and wj_. With our 
technique, fluctuations in the energy input to waves moving in opposite directions are independent. 

The aspect ratio of our simulation box is chosen to match the anisotropy of the turbulence. In 
most of our simulations, L z » L x = L y . We scale lengths to L z = 1 and velocities to va = 1- Thus 
waves take At = 1 to propagate the length of the box. The excitation level, set by the parameter 




(70) 



4.5. Simulation Design 
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(vl ± L z )/(vaL±), is chosen so that the longest waves cascade in less than At = 1. An equivalent 
statement is that a typical fieldline wanders by more than L± in the transverse direction during its 
passage across the length L z of the box. This requires the excitation parameter to be somewhat 
larger than unity. Typical values in our simulations are of order 5. We generally run our simulations 
for a few crossing times. Thus nonlinear interactions are fully expressed on all scales. 

Our basic procedure comes with a variety of refinements. The fields can be decomposed into 
their w-p and W| components as given by equations (55). Each of these may be further separated 
into shear Alfven and slow modes according to equations (56) and (57). In this manner we can 
selectively input and remove waves of any type and with any direction of propagation. 

5. SIMULATION RESULTS 

Parameters of the simulations referred to in this section are listed in Table 4. Each simulation 
is carried out in a box of dimensions L x = L y = 2 x 10~ 3 and L z = 1, includes an external 
magnetic field of unit strength aligned with the z-axis, and uses a fourth order hyperviscosity. The 
dimensionless forcing power is denoted by "P. It is chosen so that the rms dimensionless fluctuations 
of v and b have magnitude 3 x 10 -3 . These values also characterize the initial states of simulations 
of decaying turbulence. 

5.1. Power Spectra 

We obtain power spectra from our simulations as described in §4.2. 

5.1.1. ID Power Spectra 

Examples of ID power spectra obtained by averaging results from three simulations (F2, F3, 
and F4) of resolution 128 x 128 x 512 are presented in Figure 3. Each spectrum has an inertial range 
slope of approximately 1.5. The power spectra displayed in Figure 4 come from a single simulation 
(F5) with resolution 256 x 256 x 512. Aside from their extended inertial ranges, they look similar 
to those plotted in Figure 3. 

5.1.2. 2D Power Spectra 

Figure 5 displays a sequence of ID power spectra made by taking cuts parallel to the s z 
axis across the 2D power spectrum of shear Alfven waves from simulation F5. Note that there 
is negligible power at the highest s z even for the highest s±. Thus this simulation has adequate 
resolution along the z direction, something we verify for each of our simulations. As we emphasize 
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Fig. 3. — i.D Averaged Power Spectra. Energy spectra obtained by averaging results from simulations 
F2, F3, and F4 with resolution 128 x 128 x 512. 
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Alfven up wave spectrum 
Alfven down wave spectrum 
Slow up wave spectrum 
Slow down wave spectrum 
Passive scalar spectrum -10~ 7 




Fig. 4. — Highest Resolution ID Power Spectra. Energy spectra obtained from simulation F5 with 
resolution 256 x 256 x 512. 
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in §4.1 and §4.2, these cuts do not suffice to determine the structure of turbulence parallel to the 
local magnetic field. For that we need to use structure functions (see §5.2). 

5.2. Structure Functions 

Figure 6 displays transverse and longitudinal structure functions for both shear Alfven and 
slow waves calculated as described in §4.3 from data obtained by averaging results from simulations 
F2, F3, and F4. The plots are truncated at X/L = 0.5 because at greater separations the structure 
functions are affected by the application of periodic boundary conditions. 

5.2.1. Anisotropy 

Structure functions are the best measure of the scale dependent anisotropy of an MHD cascade. 
Ordered pairs of Ay and A^ obtained by equating the longitudinal and transverse structure functions 
for shear Alfven waves shown in Figure 6 are plotted in Figure 7. We leave it to the reader to judge 

2/3 

the degree to which this supports the prediction by GS that Ay oc A ± in the inertial range. 

5.2.2. Ratio of Nonlinear to Linear Time Scales 
The quantity x = (^±va)/(Mv\±) is the ratio of the nonlinear to linear timescale associated 

1 /2 

with wave packets of dimensions (Aj_, Ay). To evaluate x we take v\ ± = SFT^ (Aj_) from Figure 
6 and Aj_/Ay from Figure 7. The plot in Figure 8 establishes that x maintains a value near unity 
throughout the inertial range. 

5.3. Energy Loss 

Total, mechanical plus magnetic, energy is conserved in the inertial range of MHD turbulence. 
It is lost from high k± modes by a combination of hyperviscous dissipation and dealiasing. 13 Nei- 
ther represents reality, but we hope that their effects do not compromise inertial range dynamics. 
Figure 9 includes plots from simulation F2 of the hyperviscous and dealiasing energy losses per 
computational timestep. For reference, the total power spectrum is also displayed. Hyperviscous 
dissipation dominates dealiasing except at the highest k± where the residual power is negligible. 
This situation is typical of all our simulations. 



Energy is lost during dealiasing when we set the amplitudes of modes with \s a \ > N a /3 to zero. 
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Fig. 5. — Cuts Across 2D Power Spectrum. The Alfven energy spectrum as a function of s z at fixed 
s± from simulation F5. 




Fig. 6. — Transverse and Longitudinal Structure Functions. Structure functions transverse and 
longitudinal to the local magnetic field direction are obtained by averaging results from simulations 
F2, F3, and F4 with resolution 128 x 128 x 512. 
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Fig. 7. — Ordered Pairs of A^ and Ay. Anisotropy of MHD turbulence is quantified by plotting 
values of X± and Ay obtained by setting SFTa(X_l) = SFLa(\\) using the data displayed in Figure 
6. 
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Fig. 8. — Critical Balance. Data from Figures 6 and 7 are combined to form \, the ratio of linear 
to nonlinear time scales. Note that x has a nearly constant value close to unity throughout the 
inertial range. This confirms that MHD turbulence maintains a state of critical balance. 
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Fig. 9. — Energy Loss Per Timestep by Hyperviscous Dissipation and Dealiasing. Energy is lost 
from high k± modes by hyperviscous dissipation and dealiasing. The former dominates the latter 
when integrated over the spectrum. Neither is significant in the inertial range. The results shown 
here are from simulation F2. 



-29- 



5.4. Imbalance 

Because only oppositely directed waves interact, turbulent cascades tend to become unbal- 
anced. By unbalanced, we mean that unequal fluxes of energy propagate in opposite directions 
along the magnetic field. 

5.4-1- Forced Turbulence 

Mode energies from simulation Fl of forced turbulence with resolution 64 x 64 x 256 are plotted 
as a function of time in Figure 10. 14 Characteristic fluctuations of order unity occur on a time scale 
At = 1. Imbalance appears to saturate on longer time scales. 

5-4-2. Decaying Turbulence 

Imbalance is more severe in decaying turbulence. Figure 11 displays energies of individual 
modes as a function of time obtained from simulation Dl of resolution 64 x 64 x 256. The initial 
imbalance increases without limit. 

5.5. Passive Role Of Slow Waves 

5.5.1. Cascading Of Slow Waves By Shear Alfven Waves 

Simulation D2 of decaying turbulence with resolution 64 x 64 x 256 is designed to assess the 
mutual effects of shear Alfven waves on slow waves and vice versa. We initialize it by removing the 
upward propagating slow waves and the downward propagating shear Alfven waves from simulation 
Fl at t = 6.6. It is then run for At = 1. Figure 12 illustrates the evolution of the energies in upward 
propagating shear Alfven waves and downward propagating slow waves. The only change in the 
spectrum of shear Alfven waves is a decay at large k± which is entirely attributable to energy loss 
by hyperviscosity and dealiasing. By contrast, the spectrum of slow waves decays at all k± at a 
rate consistent with that shown in Figure 11. These findings demonstrate that shear Alfven waves 
control the MHD cascade and that the slow waves play a passive role. 



This simulation is the source of initial conditions for many higher resolution simulations. 



10 



• • • • Alfven up wave energy 

• »•• Alfven down wave energy 
Passive scalar energy 




Time 

Fig. 10. — Forced Turbulence. Mode energies as a function of time for forced turbulence from 
simulation Fl of resolution 64 x 64 x 256. 
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Fig. 11. — Decaying Turbulence. Energy as a function of time for shear Alfven and slow modes in 
decaying turbulence. The simulation is Dl with resolution 64 x 64 x 256. 
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Fig. 12. — Passive Role of Slow Waves. Upward moving shear Alfven waves interact with downward 
moving slow waves. The power spectrum of the former decays only at large k± whereas that of the 
latter decays at all k± . Results are taken from simulation D2 of decaying turbulence with resolution 
64 x 64 x 256. 
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5.5.2. Conversion Of Shear Alfven Waves To Slow Waves 

Simulation D3 of decaying turbulence with resolution 64 x 64 x 256 is tailored to measure 
the rate at which shear Alfven waves are converted into slow waves. It is initialized from Fl at 
t = 6.6 by removing all slow waves and then run for At = 1. As demonstrated by Figure 13, at 
the end of this interval, which corresponds to about a decay time at the outer scale (see Fig. 11), 
the slow waves carry negligible energy. The small admixture shown may result from the limited 
ability of our scheme of spectral decomposition to distinguish slow waves from shear Alfven waves 
as discussed in §4.1. 

5.6. Cascade Diagnostics 

We design special simulations to exploit the fact that only oppositely directed waves interact. 
Each of these is initialized by removing all but a narrow band in k± of up waves from a fully 
developed forced simulation. These are then run without forcing so that we can observe the 
evolution of the energy in the up band as it spreads into adjacent bands. Since the down waves 
evolve weakly, we restrict the lengths of these runs to At = 1/2 so that interactions do not repeat. 

Initial conditions for simulations D4, D5, D6, and D7 are provided by band-filtering simulation 
F2 at t = 2.8 with up modes retained from 2 < s± < 4, 4 < s± < 8, 8 < s± < 16, and 16 < s± < 32, 
respectively. Each of these simulations has resolution 128 x 128 x 512. Resolution 256 x 256 x 512 
simulations D8 and D9 are initialized from simulation F5 by band-filtering at t=2.95 with up modes 
retained from 16 < s± < 32 and 32 < s± < 64, respectively. 

5.6.1. Absence Of An Inverse Cascade 

Figure 14 summarizes how energy spreads from each of selected band into adjacent bands. 
It establishes that the predominant movement is toward higher k±. There is no evidence for an 
inverse cascade. A more detailed demonstration for the selected band 8 < s± < 16 is provided in 
Figure 15 which is based on simulation D6. 

5.6.2. Resolution Dependence 

Figure 16 compares results from simulation D7 of resolution 128 x 128 x 512 with those from 
simulation D8 of resolution 256 x 256 x 512. Each simulation is initialized with energy in up waves 
confined to the band 16 < s± < 32. Note how well the energies in the central and left-hand 
bands from the two simulations match as they evolve. This establishes that the largest k± band in 
simulation D7 is a valid part of the inertial range. 
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Fig. 13. — Negligible Conversion of Shear Alfven to Slow Waves. Slow wave production in a 
simulation of decaying turbulence initialized with pure shear Alfven waves. After one decay time 
the slow waves contain < 10~ 4 of the total energy. This amount is indistinguishable from the 
false slow waves that our spectral decomposition procedure would report due to the tilt of the local 
magnetic field relative to the global z axis. Data plotted comes from simulation D3 with resolution 
64 x 64 x 256. Simulation Fl provides the initial conditions for simulation D3. 
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A: 2 < s ± < 4, simulation D4 

B: 4 < s ± < 8, simulation D5 

C: 8 < s ± < 16, simulation D6 

D: 16 < s ± < 32, simulation D7 

E: 32 < s ± < 64, simulation D9 
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Fig. 14. — Summary of Bandpass Filtered Simulations. We plot the energy as a function of time 
in the center band and in the bands immediately to its left and right for each bandpass filtered 
simulation. Initially all of the energy is in the central band. As time passes it spreads into adjacent 
bands. Points correspond to simulation D8, which is initialized with the same up mode band as 
simulation D7 but with twice the transverse resolution in the down modes. There is good agreement 
between simulations D7 and D8 as shown in more detail in Figure 16. 
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Fig. 15. — Absence of an Inverse Cascade. Dashed lines depict the initial and final down wave 
spectra. The up wave spectra are plotted at a succession of times differing by At = .05. Almost 
all of the energy that leaves the band 8 < s± < 16 moves to higher s±. These data are taken from 
simulation D6. 
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Fig. 16. — Comparison of Simulations at Different Resolution. Increased resolution has little effect 
on the evolution of energy in the left-hand and central bands. Thus the latter resides in the inertial 
range in even the lower resolution simulation. Energy which moves into the right-hand band is 
more rapidly dissipated in the lower resolution simulation and more effectively stored in the higher 
resolution one. 
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5.6.3. Cascade Time 

As a standard measure of the time scale for energy transfer across Aj_, we take t c ~ Aj_/t>A x , 
where v\ ± is obtained from the transverse structure functions of downward propagating Alfven 
waves according to 2v\ = SFTai(X±). Banded simulations also permit a more direct measure of 
the cascade time as that at which the energy in the right-hand band matches that in the central 
band. We identify this version by the symbol th- 

Values for the different types of cascade time are given in Table 2. Even for the lowest k± band, 
each is substantially smaller than the time, At = I, that waves take to cross the computational 
box. 

For ease of comparison, we plot the tabulated values against k±L±/(2ir) in Figure 17. Note 
that th declines much more slowly with increasing k± than t c does. 

5.7. Intermittency 

Simulations of hydrodynamic turbulence exhibit structure that is not seen in random phase 
realizations of velocity fields with identical power spectra (Jimenez, Wray, & Saffman 1993). We 
find the same to be true for MHD turbulence. This is illustrated in Figures 18-21. The left-hand 
panels display magnitudes of the curls of upward and downward propagating shear Alfven and slow 
waves in a (x, y) slice at z = taken from simulation F5. Randomizing the phases of the Fourier 
coefficients used to generate the left-hand panels yields the images shown in the right-hand panels. 
Coherent structures, which are conspicuous in the former, are absent in the latter. 

While the eye does an excellent job of recognizing intermittency, it is helpful to have a quantita- 
tive measure. To accomplish this, we apply a sequence of high-pass filters to the Fourier coefficients 
of the Elsasser fields and a sequence of low-pass filters to their gradients. A filter is identified by 
a value of s±. High-pass filters remove modes with smaller s± and low-pass filters remove those 

Table 2. Cascade Times 
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Fig. 17. — Cascade Times. A comparison of cascade times based on different definitions. See text 
for details. 
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with larger s±. Transverse structure in the Elsasser fields is dominated by low modes and that 
in their gradients by high s± modes. Applying a sequence of high pass filters with increasing s± 
to the Fourier coefficients of the Elsasser fields emphasizes structure of decreasing scale. Similarly, 
applying a sequence of low pass filters with increasing s± to the Fourier coefficients of the gradients 
of the Elsasser fields targets structure of increasing scale. 

Filtered data is obtained from the simulations used to produce Figures 18-21. Normalized 
fourth order moments of relevant quantities q are computed according to 

M M = ^p. ( 71 ) 

where angular brackets denote volume average. 

Figure 22 displays moments of the Elsasser fields as a function k± for high-pass filtered data. 
Moments of gradients of the Elsasser fields as a function of k± for low-pass filtered data are plotted in 
Figure 23. For comparison, each figure includes moments obtained from the random phase versions 
of the corresponding simulations. It is worth noting that M^q) = 1 + 2/n for data obeying n- 
dimensional Gaussian statistics, and that slow waves correspond to n = 1 and shear Alfven waves 
to n = 2 in the limit k z . 



5.7.1. Passive Scalar 

Intermittency also characterizes the concentration of the passive scalar. In the left and right 
hand panels of Figure 24, we plot the magnitude of the gradient of the passive scalar computed in 
our highest resolution simulation F5. The contrast between the simulation and random phase data 
is striking. Coherent structures which are prominent in the former are absent from the latter. 
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Fig. 18. — Alfven Up Modes. The left-hand panel shows a grey scale image o/|V x A|| in a (x,y) 
slice at z = 0. For comparison, an image based on the same Fourier coefficients with random phases 
is shown in the right-hand panel. 




Fig. 19. — Alfven Down Modes. The left-hand panel shows a grey scale image of |V x AjJ in a 
(x,y) slice at z = 0. For comparison, an image based on the same Fourier coefficients with random 
phases is shown in the right-hand panel. 
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Fig. 20. — Slow Up Modes. The left-hand panel shows a grey scale image of | V x Sj| in a (x, y) slice 
at z = 0. For comparison, an image based on the same Fourier coefficients with random phases is 
shown in the right-hand panel. 




Fig. 21. — Slow Down Modes. The left-hand panel shows a grey scale image of |V x SjJ in a (x,y) 
slice at z = 0. For comparison, an image based on the same Fourier coefficients with random phases 
is shown in the right-hand panel. 
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Fig. 22. — Normalized 4th Order Moments for Alfven and Slow Modes. Plotted points are average 
values of moments obtained from upward and downward propagating waves. The location of the 
high-pass filter is denoted by k±_ . 
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Fig. 23. — Normalized 4th Order Moments from Gradients of Alfven and Slow Modes. Gradients are 
defined as d x q, where q is one of the Elsasser fields. Moments obtained from upward and downward 
propagating waves are averaged. The location of the low-pass filter is denoted by k±. 



Fig. 24. — Passive Scalar Gradient Magnitude. Passive scalar gradient magnitude, |Vc|, from our 
highest resolution, 256 x 256 x 512, simulation F5. The image plane is z = 0. Simulation data and 
its random phase transform are plotted in the left and right hand panels, respectively. 
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6. DISCUSSION OF RESULTS 

6.1. Comparison With GS Model 

The GSI model for the inertial range cascade of strong MHD turbulence is based on two 
assumptions. 

1. Energy transfer is local in wavenumber space. 

2. Linear and nonlinear time scales maintain near equality. 

These assumptions lead to two predictions. 

1. The ID energy spectrum E(kj_) oc kj 5 ^ 3 . 

2 /3 

2. The cascade is anisotropic with energy confined within a cone fey oc fc ± . 

Results from simulations presented in §5.1 and §5.2 agree with some aspects of the GS model and 
differ with others. An analysis of the reality and meaning of the departures from the GS scalings 
is presented in §6.5. 

6.1.1. Power Spectra and Structure Functions 

The ID power spectra displayed in Figures 3 and 4 exhibit inertial range slopes, m ps , closer 
to —3/2 than to the —5/3 predicted in GSI. This is consistent with the slopes of the transverse 
structure functions m s f shown in Figure 6 being close to 1/2. Since power spectra and structure 
functions are related by Fourier transforms, these slopes should satisfy m ps + m s f = — 1. 

A clear increase of anisotropy with decreasing scale is demonstrated in Figure 7. It is consistent 

2 /3 

with the prediction by GSI that Ay oc Aj_ . Cho & Vishniac (2000) contains the initial confirmation 
of this relation. 



6.1.2. Critical Balance 

Equality of linear and nonlinear time scales, also known as critical balance, predicts that 
Ay/fA ~ A_i_/va ± . Figure 8 shows that the ratio x = (\\ v x ± )/(^± v a) maintains a value near 
unity throughout the inertial range as predicted in GSI. However, there is a marginal problem of 
consistency. Together, Ay oc A^ 3 and \ = constant imply v\ ± oc A^ 3 . But the transverse structure 
function from which we obtain v\ ± to use in forming x yields v\ ± oc A^ 4 . 
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6.1.3. Cascade Times 



Two measures of the cascade time are plotted against k± in Figure 17; t c ~ X±/v± and 
th ~ v ±/ e - Each exhibits a power law dependence upon k± with the former's slope being steeper 
than the latter's. For v± oc A^ 3 , both t c and th would be proportional to A^ 3 . However, they 
are not. Clearly, v± oc A^ 4 , which yields t c oc A 3 / 4 and th oc X 1 / 2 , provides a better, although still 
imperfect fit. A speculative explanation for the difference between t c and th is offered in §6.6.4. 



The passive role played by slow waves in nearly transverse MHD cascades is neatly illustrated 
by Figure 12. GSII anticipates this behavior and offers a brief motivation. We provide an intuitive 
explanation in terms of field line geometry in §3.3.2. A mathematical derivation based on the 
equations of motion written in terms of Elsasser variables (eqn. [19]) is outlined below. Consider 
the evolution of upward directed waves in a cascade whose anisotropy is measured by the scale 
dependent angle ~ k^/k± <C 1. The nonlinear terms wj_ • Vwj and VP in equation (20) are 
responsible for their cascade. For comparable magnitudes of slow and shear Alfven waves, w^- Vw„ 
is smaller by a factor 6 if oc s than if oc a. Here s and a are the unit polarization vectors 
of slow and shear Alfven waves as defined in equation (18). Since the • Vw„ term is the sole 
source of P, the same comparison applies to the VP term. Note that these comparisons hold for 
both shear Alfven and slow w u waves. As they are independent of the degree of nonlinearity, they 
apply to the intermediate MHD cascade as well as to the strong one. 



Figure 13 demonstrates that the conversion of shear Alfven waves to slow waves is of negligible 
significance in MHD cascades. GSI contains the original prediction. A modified version of the 
argument given there is described below. It compares the rate at which slow waves are created in 
a balanced cascade composed entirely of shear Alfven waves to the rate at which the shear Alfven 
cascade. 

Our starting point is the Fourier transformed equation of motion for upward propagating waves 
written in terms of Elsasser variables 



6.2. 



Slow Modes 



6.2.1. Their Passive Role In Cascade 



6.2.2. Lack Of Conversion Of Shear Alfven Waves To Slow Waves 




: 2 {w T (ki) - k [k • w T (ki)] } [k • W| (k 2 )] <J(ki + k 2 - k) (72) 
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where w(k) = k z VA is the linear frequency of the shear Alfven and slow waves. The rates of change 
of the amplitudes of slow and shear Alfven waves with wavevector k in a cascade of pure shear 
Alfven waves are given by 



^-Mk))S,(k) 



8^1 



d 6 h d A k 2 



s(k)-A T (ki) ki-A^ka) <5(k!+k 2 -k), 



(73) 



and 



— - icj(k) J A T (k) 



8^/ 



d 3 h d 3 k 2 



a(k)-A T (ki) ki-A^ka) <5(ki+k 2 -k). 



(74) 



Let us compare these two rates. 15 The magnitude of s(k)-Wf(ki) is smaller than that of a(k)-w-f(ki) 
by the scale dependent anisotropy factor 6 ~ k z /k± <C 1. Thus only a fraction Q 2 (k±) <C 1 of the 
energy in shear Alfven waves is converted into slow waves as the shear Alfven waves cascade across 
k±. This accounts for the negligible production of slow waves as shown in Figure 13. 



6.3. Dynamics Of Imbalance 

The proclivity of MHD cascades for imbalance is a consequence of nonlinear interactions being 
restricted to collisions between oppositely directed waves. 



6.3.1. Forced Turbulence 



Large fluctuations are observed in the energies of different wave types in simulations of forced 
turbulence. Nevertheless, the imbalance appears to be bounded. Figure 10 provides an excellent 
example of this behavior. 

A simple dynamical model suffices to capture the essence of imbalance in forced MHD turbu- 
lence. It consists of the coupled equations 



dE^ 
~dt 



L 



1/2 



+ A, and 



dt 



L 



+ A. 



(75) 



15 The net growth rate of shear Alfven waves vanishes in a steady state cascade. Restricting the integral to fci < k 
yields the rate at which the amplitude of Af grows due to the cascading of longer (fci < k) upward propagating shear 
Alfven waves. 
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Here E denotes the energy density of the shear Alfven waves, L the transverse outer scale, and 
A the excitation rate. The equilibrium energy density and the nonlinear cascade time scale are 
defined by E eq = (AL) 2 / 3 and t c = L^/A 2 / 3 . 

To investigate the stability of forced balanced cascade, we set 

£ T =£ cq + A£ T and E± = E^ + AE±, (76) 

and then substitute these expressions into equations (75) to obtain 

^ = -i-(2A £t + A £l ) and ^ - --L (2AE, + AE,) . (77) 

Assuming a time dependence proportional to e st , we find eigenvalues 

si = ~ and s 2 = -±-. (78) 

This establishes the stability of the forced balanced cascade. It also shows that fluctuations asso- 
ciated with the si eigenmode decay rather slowly. These characteristics accord well with the runs 
of Alfven wave energy densities displayed in Figure 10. 

A more sophisticated analysis would include a proper statistical treatment of forcing and an 
investigation of the spectrum of fluctuations. 



6.3.2. Decaying Turbulence 

Simulations of decaying MHD turbulence exhibit large imbalances. Thus a perturbation anal- 
ysis is inappropriate. Fortunately, for A = equations (75) admit an analytic solution. As is easy 
to verify by direct substitution, AE 1 / 2 = E 1 ^ 2 — E 1 ^ 2 = AE^ 2 is a constant, and that 

Thus imbalance grows exponentially in decaying turbulence. This accounts qualitatively for the 
behavior seen in Figure 11. 

Dobrowolny, Mangeney, & Veltri (1980) propose the growth of imbalance in decaying MHD 
turbulence as an explanation for the fact that the preponderance of shear Alfven waves in the 
solar wind propagate outward along the interplanetary magnetic field. Support for this proposal is 
provided by simulations described in (Pouguet, Meneguzzi, & Frisch 1986). 



6.3.3. Axial Asymmetry 



Axial asymmetry refers to a net polarization of Shear Alfven waves. This can occur even in a 
balanced cascade. MHD cascades have a tendency to develop axial asymmetry because the strength 
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of nonlinear interaction between oppositely directed shear Alfven waves 1 & 2 is proportional to 
a(ki) -a(k2) (see eqn. [74]). Thus the interaction vanishes for parallel polarizations and is strongest 
for orthogonal polarizations. 

Decaying turbulence is unstable to the growth of axial asymmetry. Waves with the subdomi- 
nant polarization cascade more rapidly that those with the dominate polarization. Axial asymmetry 
is bounded in forced turbulence. However, correlated fluctuations occur across the inertial range 
within regions of spatial scale comparable to the outer scale since the cascade tends to preserve 
polarization alignment. 

6.4. Intermittency 

Tubes of high vorticity, often referred to as worms, are prominent features in simulations 
of hydrodynamic turbulence. Worms have diameters of order the dissipation scale and lengths 
approaching the outer scale. They are thought to form from the rolling up of vortex sheets. In 
spite of their prominence, worms do not affect the inertial range dynamics (Jimenez, Wray, & 
Saffman 1993). 

Coherent structures are also evident in MHD simulations. Examples are shown in §5.7 where 
the magnitudes of the curls of the dynamical fields and of the gradient of the passive scalar are 
plotted in (x, y) slices. These regions have narrow dimensions comparable to the dissipation scale 
and lengths approaching the outer scale L±. In these respects they resemble worms. We suspect 
that these structures are vortex sheets which extend along the z axis and that the magnetic field 
prevents them from rolling up. Their correlation lengths along z are unresolved. However, this is 
not a strong constraint since kM z Lj_ rs 0.3. 

6.5. Comparison With Previous Simulations 

Shebalin, Matthaeus, & Montgomery (1983) report the development of anisotropy in isotropi- 
cally excited MHD. Their simulations are two-dimensional, with one axis parallel to the direction 
of the mean magnetic field. Thus they are composed entirely of slow waves. An isotropic distri- 
bution of slow waves will initiate an anisotropic cascade. However, nonlinear interactions among 
slow waves weaken as the cascade becomes more transverse because their strength is proportional 
to coefficients such as k2 • S^(ki), and for k — ► kj_, s — > z. 16 Convincing demonstrations of the 
development of anisotropy in fully three-dimensional MHD simulations are presented in Oughton, 
Priest, & Matthaeus (1994) and in Matthaeus et al. (1998). Each of these papers provides evi- 
dence that anisotropy increases at smaller scales. Each also claims that up to a saturation limit, 



This remains true in 3D 
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anisotropy is more pronounced the larger the ratio of mean to fluctuating magnetic field strength. 
Matthaeus et al. note that this latter trend is inconsistent with the scaling for anisotropy proposed 
by GSI. 

A analysis by Cho & Vishniac (2000) clears up the confusion regarding the scale dependence 
of anisotropy in MHD turbulence. Unlike previous workers, who measure anisotropy in coordinate 
systems fixed to the sides of their computational boxes, Cho and Vishniac compute anisotropy in 
local coordinate frames tied to the direction of the total magnetic field. We follow their technique 
of using structure functions computed in directions parallel and perpendicular to that of the local 
magnetic field (see §5.2.1). Both they and we find results that are consistent with the relation 

2/3 

Ay oc A ± proposed by GSI. However, the ratio of the fluctuating to mean field is ~ 0.5 in their 
simulations and ~ 0.01 in ours. 



6.6. Departure From GS Scalings 

The GS scalings lead to the unambiguous prediction of a scale dependent anisotropy, ku oc , 

—5/3 

and a ID Kolmogorov spectrum, E(k±) cx kj_ . While our simulations are in accord with the 
former, they consistently indicate that the ID power spectrum has an index closer to 3/2 than 5/3. 

We do not know whether our simulation method is producing an anomalous power law or 
whether we have discovered a feature of MHD turbulence that is not incorporated in the GS 
scalings. We first examine several effects which might produce anomalous power laws in numerical 
simulations. Then we offer some speculations about the role of intermittency. 



6.6.1. Forcing 

Borue & Orsag (1994) attribute the fc _L85 inertial range they find for a simulation of hydro- 
dynamic turbulence to temporal intermittency associated with forcing. We are unable to discern 
any difference in the inertial range slopes from our simulations of forced and decaying MHD tur- 
bulence. As an example, compare the shear Alfven wave spectra from both forced simulation Fl 
and decaying simulation D3 which are plotted in Figure 13. 



6.6.2. Dealiasing 

Energy loss in our simulations is due to a combination of dealiasing and hyperviscosity as illus- 
trated in Figure 9 for F2, one of our intermediate resolution simulations. Their relative importance 
depends upon the parameter V(ttN/2) 1 ^ 3 defined in equation (68). With our choice of viscosity pa- 
rameters, dealiasing contributes relatively less of the energy loss in our lowest resolution simulation 
Fl and relatively more in our highest resolution simulation F5. However, all of our simulations 



-52- 



show similar inertial range slopes as a comparision of Figures 3, 4, and 13 verifies. If there is any 
significant difference, it is that F5, the highest resolution simulation, for which dealiasing is least 
important, has the steepest inertial range slope. 

6. 6. 3. Hyperviscosity 

Hyperviscosity applied in hydrodynamic simulations is known to produce a spurious flattening 
of the ID inertial range slope over a range of k approaching the viscous cutoff (Borue & Orsag 
1994). A weaker form of this bottleneck effect is apparent in the simulations of MHD turbulence 
described by Cho & Vishniac (2000). Both Borue k Orsag (1994) and Cho & Vishniac (2000) 
employ an 8'th order hyperviscosity. Our simulations, which use a 4'th order hyperviscosity, show 
no indication of a bottleneck effect. This can be seen from the absence of flattening of the inertial 
range slope at high k± in Figures 3 and 4. 

Could hyperviscosity flatten the slope across the entire inertial range? Muller & Biskamp 
(2000) suggest that it does. They present the results of a 512 3 simulation with nearly isotropic 
forcing which uses l'st order (physical) viscosity. It exhibits an inertial range slope slightly steeper 
than 5/3. Then they mention that a similar calculation done with 2'nd order hyperviscosity results 
in a flatter inertial range. To test the effect of hyperviscosity on the inertial range slope, we carry 
out simulation F6 which uses physical viscosity but is otherwise similar to simulation F2. The ID 
spectrum from this simulation is plotted in Figure 25. Its inertial range appears to have a slope 
closer to 3/2 than 5/3. This is not entirely conclusive because the inertial range is truncated at 
the low k±_ end by forcing and at the high k± end by viscosity. Higher resolution simulations which 
are beyond our current computational resources are needed to settle this issue. 

6.6.4- Speculations 

For the moment, let us accept the shallow inertial range slope as a real feature of MHD turbu- 
lence. How might it be accounted for? An intriguing possibility is that the nonlinear interactions 
responsible for the cascade become increasingly intermittent with decreasing scale. 

A given degree of spatial intermittency in the energy density is likely to have more serious 
consequences for the turbulent cascade in MHD than in HD. In HD energy cascades according 
to the local value of X/v. But in MHD nonlinear interactions are restricted to collisions between 
oppositely directed wave packets. Thus if the spatial filling factor of the energy density, /, is small, 
that of the interactions, / 2 , is smaller still. This may account for the shallower slope of t^ ~ v"j_/e 
as compared to t c ~ \±/v± seen in Figure 17 and discussed in §6.1.3. 

Research reported in this paper was supported by NSF grant 94-14232 and a NSF Graduate 
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k ± L ± / (2tt) 

Fig. 25. — Inertial Range Computed With Ordinary Viscosity. Simulation F6 with v\ = 2 x 1CP 4 
and resolution 128 x 128 x 512 is similar to simulation F2 but uses ordinary viscosity instead of 
hyperviscosity. 
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A. APPENDIX 

A.l. Spectral Notation 

Cartesian coordinates are distinguished by Greek indices which run from 1 — 3. Simulations 
are carried out in boxes whose sides have lengths L a which are partitioned by N a grid points. 
Integer coordinate components, l a , and integer wavevector components, s a , are defined through 
the relations 

x a = -rj-la, where < l a < N a , (Al) 

and 

k a = —s a , where - — < s a < — . (A2) 

L/ a Z Z 

The discrete Fourier transform in ID is given by 

i 

where the tilde, ~ denotes Fourier transform. Generalization to 3D is trivial. 



A. 2. Spectral Algorithm 

A. 2.1. Fourier Space Equations 

We evolve the incompressible MHD equations in Fourier space where they take the form 
(Lesieur 1990) 

d t v a = -iky (^5 a(3 - ^) {vpu^ - bpb^j - u n k 2n v a , (A4) 

d t b a = -ik l3 (v f3 b a - bpvj - u n k 2n b a , (A5) 

k a v a = 0, (A6) 

k a b a = 0, (A7) 

d t c = —ikpvpc — v n k 2n c. (A8) 
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A.2.2. Integration Method 

Equations (A4), (A5), and (A8) constitute a system of ordinary differential equations with 
time as the dependent variable and the Fourier coefficients {v a , b a , 5} as the independent variables. 
We employ a modified version of the second order Runge-Kutta algorithm (RK2) to advance the 
variables in time. First order algorithms are substantially less stable than RK2 at the same timestep. 

RK2 advances the variables across an interval At in two stages. Derivatives evaluated at the 
initial time are used to compute trial values of the variables at the midpoint At/2. Then derivatives 
computed at At/ 2 with these trial values are used to advance the variables from t = to At. In 
symbolic form 

g tria i (At/2) = g(0) + dtq(0) At/2 (A9) 

is followed by 

q(At) = q{0) + d t q{At/2)At, (A10) 

where dtq(At/2)At is evaluated using q tr \ a \(At/2). Each stage involves a first order Euler (El) step 
in which the derivative is taken to be constant. 

We make one departure from standard RK2 and treat diffusive terms with an integrating 
factor. Consider an equation of the form 

d t q(k)=A-v n k 2n q(k), (All) 



where A comprises the non-diffusive terms. Its solution, with A constant throughout the interval 
At, is 

-v n k 2 "At ( M2 ) 



q(At) 



m + ^(^ M - 1) 



We use this expression in place of El in each stage of RK2. To lowest order in v n k 2n At, equation 
(A12) reduces to El. However, it has the advantage that it yields stable solutions to equation 
(All) with constant A for arbitrary values of v n k 2n At whereas El yields unstable solutions for 
v n k 2n At > 2. 



A. 2. 3. Dealiasing 

Bilinear terms in equations (A4), (A5), and (A8) are calculated by transforming the individual 
fields to real space, carrying out the appropriate multiplications there, and then transforming the 
products back to Fourier space. This requires N1N2N3 log Ni log N2 log N3 operations using the 
Fast Fourier Transform (FFT) algorithm; (N1N2N3) 2 operations would be needed to carry out the 
equivalent convolution in Fourier space. 

This economy comes at the price of either a 1/3 reduction in resolution or an aliasing error 
(Canuto 1988). To appreciate this, consider the ID product 

PQ(s) [£✓ p(s')e- 2ms ' l / N J2 S „ q{s")e- 2ms " l / N , ! > J A 



-56- 



= J? E s > P^)q^')^ + '" )l/N S 8 ,s'+s''+ m N, (A13) 

where m is any integer. The m = terms comprise the convolution, and the remainder the aliasing 
error. To avoid the aliasing error, we set all Fourier components with |s| > N/3 to zero both before 
we compute the real space fields and again after we return the bilinear terms to Fourier space. 
Truncation ensures that Fourier components of bilinear terms with m ^ vanish. Its cost is the 
reduction of the effective spatial resolution from N to 2N/3. 



A.3. Tests Of The Spectral Code 

Time derivatives of field quantities computed with the spectral code agree with those obtained 
from a finite difference program with an elliptic incompressible pressure operator. Although the 
latter is unstable, it offers an independent method for computing time derivatives. 

The code preserves the solinoidal character of v and b. To machine accuracy it returns 
k • cVv(k) = and k • 9jb(k) = 0. It also conserves energy. Provided v n = 0, it yields dtE = 
^ k v(k) • <9(v(k) + b(k) • d t b(k) = 0, again to machine accuracy. 

Harmonic Alfven waves are evolved by our spectral code in a manner consistent with their 
analytic dispersion relation. 

Results obtained from a simulation of decaying hydrodynamic turbulence (Z45) run with our 
code agree with those from a more thorough simulation by Jimenez, Wray, & Saffman (1993). Our 
simulation is carried out in a cubic box with L = 1.0, has resolution 256 3 , kinematic viscosity 
v = 8 • 1CT 4 , timestep At = 2.5 x 1CT 4 , and is initialized with rms velocity v = 1.0. We compute 
components of the velocity gradient longitudinal, Vy^y, and transverse, Vj_i>||, to v at each point 
in our computational box. Distribution functions, PF q (x), of each quantity, q, are compiled and 
moments calculated according to 



"OO 

M n = . J -°°~ " g 7 B/2 - (A14) 



These are shown in Table 3. 



Because our simulation is of decaying turbulence whereas that of Jiminez et al. is forced, appropriate 
comparisons are restricted to inner scale quantities derived from components of Vv and exclude 
outer scale quantities derived from components v. With this proviso, our results are in satisfactory 
agreement with theirs. 



Table 3. Tests of Spectral Code 



Simulation Z45 Jiminez et al. Gaussian 



n 


V 


V||f|| 


V±v\\ 


V 


V||W|| 


VjV|| 




3 





-.43 








-.50 








4 


1.6 


4.4 


5.9 


2.8 


4.6 


6.19 


3 


5 





-6.5 








-8.0 








6 


3.4 


48 


102 


13.0 


55 


110 


15 



Table 4. Simulation Parameters 



ID 


(N ± ,N Z ) 




At 




1/4 




Comments 


Fl 


64, 256 


4 


x 10~ 4 


5 


x 


10~ 


-37 


V = 2 x 10~ 5 


F2 


128,512 


3 


x 10- 4 


5 


x 


10~ 


-40 


V = 2x 10~ 5 


F3 


128,512 


3 


x 10- 4 


5 


x 


10~ 


-43 


V = 2x 10~ 5 


F4 


128,512 


3 


x 10~ 4 


5 


X 


io- 


-43 


V = 2x 10~ 5 


F5 


256,512 


1. 


5 x 10~ 4 


5 


X 


io- 


-43 


V = 2x 10~ 5 


Dl 


64, 256 


4 


x 10~ 4 


5 


X 


io- 


-40 


V = 


D2 


64, 256 


4 


x 10~ 4 


5 


X 


io- 


-37 


V = 0, A T +S i 


D3 


64, 256 


4 


x 10- 4 


5 


X 


io- 


-37 


V = 0, A T +A [ 


D4 


128,512 


3 


x 10- 4 


5 


X 


io- 


-40 


2 < s± < 4 


D5 


128,512 


3 


x 10" 4 


5 


X 


io- 


-40 


4 < s± < 8 


D6 


128,512 


3 


x 10" 4 


5 


X 


10" 


-40 


8 < s± < 16 


D7 


128,512 


3 


x 10" 4 


5 


X 


io- 


-40 


16 < s ± < 32 


D8 


256,512 


1. 


5 x 10~ 4 


5 


X 


10" 


-43 


16 < s ± < 32 


D9 


256,512 


1. 


5 x 10~ 4 


5 


X 


10" 


-43 


32 < s± < 64 
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A.4. Catalog Of Simulations 

A. 4-1. Simulations of Forced Turbulence 

Our basic simulations include forcing at a total power per unit mass of V = 2 x 1CP 5 . Recall 
that p = 1 and va = 1- Statistically equal power is input into shear Alfven and slow waves 
propagating in opposite directions along the magnetic field. Higher resolution simulations run for 
shorter times are initialized by refining lower resolution simulations run for longer times. 

Our sequence of forced simulations begins with Fl which has resolution 64 x 64 x 256 and 
runs up to t = 6.6. Initial condition for simulations F2, F3, F4 with resolutions 128 x 128 x 512 
are drawn from Fl at t = 2.4, 4.7, and 6.6, respectively. These are times at which the fluxes of 
oppositely directed shear Alfven waves in Fl nearly balance. The refined simulations run for an 
additional At = 0.4, long enough for small scale structure to develop up to the dealiasing cutoff. 
Our highest resolution, 256 x 256 x 512, simulation F5 is initialized from F2 at t = 2.8 and run 
until t = 2.95. 

A. 4-2. Simulations of Decaying Turbulence 

Our simulations of decaying turbulence are designed to test specific properties of the MHD 
cascades. Simulation Dl continues Fl without forcing from t=2.8 to t=9.9. Simulations D2 and 
D3 are initialized from Fl at t = 6.6, the former by removing the Alfven down and slow up 
waves and the latter by removing all slow waves. A series of simulations are initialized from forced 
simulations by removing all upward propagating waves outside a specified band while leaving the 
down modes unchanged. Simulations D4, D5, D6, and D7 are each initialized from F2 at t=2.8 
with the up modes band-filtered from 2 < s± < 4, 4 < s± < 8, 8 < s± < 16, and 16 < s± < 32, 
respectively. Likewise, simulations D8 and D9 are initialized from F5 at t=2.95 with up modes 
band-filtered from 16 < s± < 32 and 32 < s± < 64, respectively. The former has the same up band 
as D7 but twice the transverse resolution. 
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